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Abstract The non-dissipative relativistic force-free condition should be 
a good approximation to describe the electromagnetic field in much of the 
pulsar magnetosphere, but we may plausibly expect it to break down in 
singular domains. Self-consistent magnetospheric solutions are found with 
field lines closing both at and within the light-cylinder. In general, the 
detailed properties of the solutions may be affected critically by the physics 
determining the appropriate choice of equatorial boundary condition beyond 
the light-cylinder. 

1 Introduction 

With canonical values inserted for the neutron star dipole moment and the 
particle density, and allowing for realistic pair production efficiency near 
the star in the wind zone, the electromagnetic energy density in the pulsar 
magnetosphere will be much greater than the kinetic energy density except 
for very high values of the particle 7- values. This suggests that over much of 
the magnetosphere, the relativistic force-free equation will be a good approx- 
imation for determining the magnetic field structure. Equally, experience 
with the analogous non-relativistic equation suggests that there are likely 
to be local domains in which non-electromagnetic forces (including possi- 
bly inertial forces) are required for force balance. An associated question 
is whether the dissipation-free conditions can be imposed everywhere, or 
whether the dynamics of the problem will itself demand local breakdown in 
the simple plasma condition E • B = in singular domains, in addition to 
the acceleration/pair production domain. 

Of the various discussions of the axisymmetric pulsar magnetosphere 
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in the literature, which go back to the early days of pulsar theory (e.g. 
Scharlemann and Wagoner 1973), we refer particularly to the recent paper 
by Contopoulos et al. 1999, from now on referred to as CKF, and to those by 
Michel 1973a,b, 1991, Mestel and Wang (MW) 1979, Mestel et al. (MRWW) 
1985, Fitzpatrick and Mestel (FM) 1988a,b , Mestel and Pryce (MP) 1992, 
and Mestel and Shibata (MS) 1994. (Much of the relevant material appears 
in Mestel 1999, 2003.) Most models have within the light-cylinder ('1-c') a 
'dead zone,' with field lines that close within the 1-c, and without any current 
flow along the (purely poloidal) field, and a 'wind zone,' with poloidal field 
lines that cross the 1-c, and with poloidal currents maintaining a toroidal 
field component. The spin-down of the star occurs through the action of the 
magnetic torques associated with the circulatory flow of current across the 
light-cylinder. 

Although we are here discussing an idealised mathematical model of 
the pulsar magnetosphere - not even allowing the pulsar to pulse! - its 
structure may nevertheless have great relevance to the class of real pulsars 
whose magnetic axis is not far from alignment. A striking example of this is 
PSR0826-34 (Biggs et al 1985), a pulsar which emits throughout the entire 
pulse period, suggesting our line of sight is almost along the rotational axis. 
Other bright pulsars which are suspected on the basis of their profile widths 
and relatively slow spin-down rates to be near alignment include PSR0031- 
07, PSR0943+10 and PSR0809+74 (see Lyne and Manchester 1988 and the 
catalogue in Rankin (1993)). 

All these pulsars exhibit highly regular subpulse 'drift', by which the 
sub-pulses march gradually and steadily through the pulse window at a 
fixed rate. This suggests that their magnetospheres are in a quasi-steady 
state, and that steady-state models such as those discussed here may have 
direct physical relevance. Furthermore, the fact that all the above pulsars 
occasionally - yet apparently spontaneously - switch to one or more alter- 
nate 'drift' rates, accompanied by a widening (or narrowing) of the pulsar 
emission profile (Wright and Fowler 1981, Vivekanand and Joshi 1997, van 
Leeuwen et al 2002) would indicate that the new steady emission has moved 
to outer (inner) field lines. Thus it would appear that more than one steady 
state is possible in the same pulsar - presumably corresponding to different 
polar cap shapes, and hence to different 'dead zone' boundaries. This would 
add relevance and motivation to our discussion in Section 6, where we derive 
magnetosphere solutions for a range of assumptions about the location and 
properties of the dead zone boundary. 
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2 The relativistic force-free equation 



The outer crust of the neutron star of radius R rotates with angular velocity 
a about the axis defined by the unit vector k. We confine attention in this 
paper to steady, axisymmetric states, and with the magnetic axis parallel 
rather than anti-parallel to the rotation axis k. The notation is as in the 
cited papers (e.g. FM). The cylindrical polar coordinate system (w,<f>,z) is 
based on k; t is the unit azimuthal, toroidal vector. The poloidal magnetic 
field Bp is everywhere described by the flux function P(w, z): 

„ n ( t \ 1 fdP dP\ 
Bp = -VP x (-)=-(_,(,,--). (1) 

By Ampere's law, B p is maintained by the toroidal current density 

j t = (c/4f)(VxB p ) = (c/4f)(VxB)^t (2) 
= {c/4tt)[V 2 P/w - (2/w 2 )dP/dw]t. (3) 

The basic, unperturbed stellar field is assumed to have a dipolar angular 
dependence and polar field strength B s . If as assumed, aR/c <C 1, it will 
be found that the magnetospheric currents have little effect on the form of 
P near the star, so that it is then well approximated locally by the vacuum 
form 

P = _^! ^ ( 4 ) 

In the wind zone of an active magnetosphere there will be also a toroidal 
component B t , conveniently written 

B 4 ,t = -[4TrS/c]{t/w), (5) 

where S clearly must vanish on the axis. By Ampere's law, the field (5) is 
maintained by the poloidal current density 

j p = ( c /4vr)V x B t = -VS x t/zu. (6) 

Thus S(w, z) is a Stokes stream function, constant on the poloidal current 
lines, with —2ttS measuring the total outflow of charge between the axis 
and the current line S. In a steady state the total outflow of charge must 
be zero, so there must be a current closing streamline on which S vanishes. 
The crust is taken to be a classical perfect conductor, with the electric 

field 

E = —{aw/c)t x B. (7) 
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In the surrounding magnetosphere, E is in general written as the sum of a 
'corotationaP part of the form (7) and a 'non-corotationaP part V^: 



E = —{aw/c)t x B — Vip, (8) 

where the function ip is determined by the magnetospheric physics. 

A systematic treatment of the dynamics of a 'cold,' dissipation- free 
electron-positron gas has been given by Melatos and Melrose (1996) (see 
also Blackman and Field 1993). The significant differences from a normal 
plasma are due in part to the particles being relativistic, and also through 
both species having the same rest-mass, instead of differing by a factor 
Amu/rric in standard notation. The two-fluid, collision-free and so non- 
dissipative equations, written in terms of 'lab-frame' number densities 
and velocities v^, and Lorentz factors 7 ± = {1 — (i^/c) 2 } -1 / 2 , are reduced 
to one-fluid equations in terms of 

n=(n+ + n-), U = + — , T = (1 - U 2 /c 2 )^ 2 

p c = e{n + — n~), j = e(n + v + — re~v~) (9) 

(The lab-frame densities n ± are related to the proper densities by n^ 1 = 
7 ± n ± .) 

The particle continuity and current continuity equations retain their 
standard forms 

^ + V-(nU) = 0; ^ + V-j = 0. (10) 
The equation of bulk motion with velocity U reduces to 

mn^l + mnV • V(TU) = ^~P^ xB (11) 
ot c 

provided j <C nell, even if, as expected in a relativistic problem, the 'con- 
vection current' p c ll due to bulk motion of the net charge density is not 
small compared with the current due to relative motion of the two species. 
The generalized Ohm's law reduces to 



m d /rj\ m TT _ /Ti\ m .. __. _,„ TN ^ UxB 



As stated, in this paper we are interested primarily in steady states that 
are force-free over the bulk of the domain. Just as in a classical plasma, in 
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an zte-plasma the total particle number density will normally greatly exceed 
the GJ number density, defined by (16) below, so the 'perfect conductivity' 
and the 'force-free' conditions are distinct: a gas satisfying the 'perfect 
conductivity' condition need not be force-free. However, one can easily check 
that if inertial terms on the left of (11) are small compared with either term 
on the right, then the inertial corrections in (12) are a fortiori small. Thus in 
the absence of collisions, radiation damping or electron-positron annihilation 
- all neglected in (11) and (12) - the force-free condition implies also the 
ideal MHD (perfect conductivity) condition. To sum up: for our problem, 
the appropriate approximations to the dynamical equations are the ideal 
MHD condition 

E + UxB/c = 0, (13) 

implying EB = 0, and the mutual cancellation of the electric and magnetic 
contributions to the Lorentz force, i.e. the vanishing of the term on the right 
of (11), 

(j - p c U) x B jxB 

^ — ^—i- = p c E + = 0, 14 

c c 

on use of (13). 

It should be noted that in the more general, non-axisymmetric case, 
with k and p non-aligned, the time-dependent terms will not be negligible 
everywhere. As discussed by Melatos and Melrose (1996), at a finite distance 
far beyond the 1-c, the displacement current will dominate over the particle 
current, and the MHD approximations must break down. 

We limit discussion explicitly to the case in which the simple plasma 
condition E • B = is supposed to hold throughout the magnetosphere, all 
the way from the rigidly rotating, perfectly conducting star out to the 1-c 
(the 'inner domain'), and beyond into the 'outer domain' between the 1-c and 
infinity. Then from (8), B • VY> = 0, and the constant value of ip through- 
out the stellar crust, implied by (7), is propagated into the magnetosphere, 
yielding 

E = -{aw/c)t x B = (azu/c)(-B z ,0,B x ) = (a/c)VP, (15) 

so that all field lines corotate with the star. This assumes that within the 
dead zone, there is no vacuum gap separating the negatively and positively 
charged domains (Holloway and Pryce 1981; MRWW; FM). In the open field 
line domain, near the star there has to be a locally non-trivial component 
of E along B p , able to accelerate the primary electrons to 7- values high 
enough for pair production to occur. If, as assumed in this paper, the 
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electron-positron plasma does achieve a steady state (cf. Shibata et al 1998, 
2002), the electric field will again satisfy E • B = but will now be given 
by E = — a{P){w / c)t x B - i.e. with the field lines beyond the acceleration 
domain having individual rotation rates a{P) that differ somewhat from the 
rotation a of the star (cf. MS, Section 4). However, at least for the more 
rapid rotators this effect will be fairly small, and will for the moment be 
ignored. 

By the Poisson-Maxwell equation, the Goldreich- Julian (1969) (GJ) 
charge density maintaining the electric field (7) is 



V-E 

47T 



2vrc 
a 

2vrc 



B--rx(VxB) 



B z - -w{V x B) 



(16) 



With the rotation and magnetic axes aligned, the primary outflowing par- 
ticles are the negatively charged electrons, yielding a negative current j p 
emanating from the polar cap. The stream function S defined by (6) be- 
gins by increasing from zero on the axis, so that is negative - the field 
lines are twisted backwards with respect to the axis k. As the electric force 
density 

p c E = -p c (aw/c)t x Bp (17) 

is purely poloidal, in a force-free magnetosphere the toroidal component of 
the magnetic force density j p x B p /c vanishes (the 'torque- free' condition), 
so j p must be parallel to B p , yielding from (1) and (6) the functional relation 

S = S(P), Jp = ^B p : (18) 

the poloidal current streamlines are identical with the poloidal field lines. 
The respective contributions of B p and B t to the poloidal force density are 



j t x B p /c = (V x B)^(t x B p )/4vr 



and 



„ , AttS dS , _ . 
j p xB t /c=- I - — (txBp) 



47T 



c 2 w 2 ~ dP 



S^VP, 



(19) 



(20) 



c 2 wdP K 

on use of (18) and (1). By (1), (16), (17), (19) and (20), the poloidal 
component of the force-free equation is 



— (V x B) 

w 



1 



aw 
c 



+ 



aw 
c 



2B Z 



w^ 



+ 



cw) dP 



0. (21) 
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Note that the first term in (21) combines part of the electric force with the 
force due to B p ; the second again comes from the electric force, while the 
third is due to B t . 

With the sign convention in (1), P decreases from zero on the axis, 
so that initially dS/dP is negative, and the force density (20) due to the 
toroidal field B t acts toward the axis. Since S must vanish on the current 
closing streamline, there will in general be an intermediate field line on which 
dS/dP and so also the volume current density j p changes sign. Below this 
field line the force (20) acts towards the equator. Note that the return of S 
to zero need not be continuous: closure via a sheet current is not excluded 
(cf. Section 4). In fact at this stage, one cannot rule out that S increases 
monotonically from zero to its maximum on a limiting field-streamline and 
then drops to zero, so that all the return current is in the sheet. The actual 
allowed form for the function S(P) will emerge as part of the global solution. 

As in the earlier work, we define dimensionless coordinates (x,z) = 
(a/c)(w, z), and normalize P, S, E,B in terms of a standard light-cylinder 
field strength B lc = (B s /2)(aR/c) 3 : 

P = PB lc (c/a) 2 , S = SB lc (c 2 /47ra), (B, E) = £ lc (E, B). (22) 

(Once defined, the dimensionless quantities are again immediately written 
without the bars). The normalized fields have the form 

B p = -VPxt/x, B 4> = -S/x, E = VP, (23) 

with P satisfying 

, 2 ,d 2 P (l + x 2 )dP , 2 x d 2 P AS , . 

As (w, z) — > 0, P must reduce to the normalized point dipolar form — x 2 / (x 2 + 
z 2 ) 3 / 2 . The light-cylinder x = 1 is a singularity of the differential equation 
for P. 

An illustrative solution of (24) within the 1-c for the non-active case 
S = has been discussed by Michel (1973, 1991), Mestel & Wang (1979) 
and Mestel & Pryce (1992), and will be referred to as the MMWP-field. The 
field lines which reach the equator cross normally, forming a closed domain: 
the inner domain equatorial boundary condition is 

B x (x,0)=0, (25) 
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with no equatorial current sheet. The flux function P is most easily con- 
structed as a Fourier cosine integral (cf the Appendix). If neither the field 
nor the volume current density are to be non-singular at the 1-c, (24) with 
S = requires that B z = —dP/dx = when x = 1: the field lines cross the 
1-c normally, and there is a neutral point at the intersection (1,0) of the 1-c 
and the equator. 

The MMWP model is of pedagogic interest through its showing how 
in this essentially relativistic problem, the macroscopic electric force acting 
on the corotating GJ charge density causes a marked deviation of B from 
the curl-free dipolar form as the 1-c is approached. However, the solution 
is clearly incomplete, as there is no treatment of the domain beyond the 
1-c. A viable, strictly inactive model will in fact have large vacuum gaps 
within the 1-c (e.g. Smith et al. 2001 and references within). As discussed 
in Sections 5 and 6, some features of the MMWP field should persist in a 
realistic active model, such as the large dead zone, and a similar departure 
from the curl-free structure near and beyond the 1-c. 

3 The domain beyond the 1-c 

Equation (13) and (15) combine into 

(U - ami) x B = 0, (26) 

yielding 

U = kB + awt (27) 

- the sum of corotation with the star plus flow parallel to the total field 
B = Bp + B t , a result familiar from standard stellar wind theory. If U p 
were to vanish beyond the 1-c, then from (9), n + (vf — am)+n~(vj~ — am) = 
0, requiring that either the electrons or the positrons would have to be 
superluminal. Thus along field lines that cross the 1-c, there is not only a 
net current but also a net particle flow, as is indeed implied by the term 
'wind zone.' Then if, as assumed, the simple 'perfect conductivity' condition 
(26) continues to hold everywhere beyond the 1-c, so that there is strictly 
no trans-field motion of the plasma, then in a steady state the field must be 
topologically 'open,' with no field lines crossing the equator. 

A CKF-type model has the field lines crossing the equator normally 
within the 1-c, as in the MMWP and FM fields, but beyond the 1-c, the 
boundary condition B z = is imposed at the equator. In the northern 
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hemisphere B x > 0, so by (7), near the equator the electric field xB x is in 
the positive z-direction. As the basic field is dipolar, B x must change sign 
at the equator, implying a locally positive current density and a magnetic 
force density that acts towards the equator. Thus the poloidal magnetic field 
pinches, as in familiar non-relativistic problems, but the combined electric 
force (17) and the B p force (19) is [(x 2 — 1)(V x B)0/47r]k and so acts away 
from the equator, since x > 1. The proposed force- free equilibrium near the 
equator must therefore be maintained by the pinching effect of B t . We have 
seen above that SdS/dP does indeed become positive at low latitudes, so 
that (20) has the required direction. 

Consider first the idealized case, with B x , B^ and E z abruptly reversing 
sign at the equator, and so with positive sheet currents Js, J x and positive 
surface charges a. Over the surface z = +e, there is a net electromagnetic 



stress Tijkj where 



T — 



13 4vT 



-- (E 2 + B 2 ) Sij + {EiEj + BiBj) 



(28) 



is the Maxwell stress tensor. The electric terms yield (x 2 B 2 /8ir)\i and the 
magnetic terms — [(B x + B^) /8ir]k, so the net electromagnetic pressure, act- 
ing towards the equator, is 

[Bl - {x 2 - l)B 2 x ]/8n. (29) 

There is an equal electromagnetic pressure, also acting towards the equator 
on the surface z = —e. Provided 



\B^,/B X \ > Vx 2 - 1, (30) 

then an equal particle pressure at the equator is both necessary and sufficient 
to maintain equilibrium. 

The force-free equation (21) can be written succinctly as dTij/dxj = 0. 
Near the equator, \B Z \ << \B X \ in this model, and the force-free condition 
becomes 



8x3 33 8tt dz 



B 2 + (x 2 -l)B 2 =0, (31) 



whence 

j-[B 2 -(x 2 -l)B 2 ]=p eq (32) 

where p eq is independent of z. If (30) holds, then p eq is identified as a particle 
pressure at the equator, where B x and Bs change sign. 



9 



More realistically, the transition to zero field at the equator will be 
continuous, with the particle pressure p increasing steadily across half of the 
thin equatorial sheet. The force- free condition (31) is replaced by 



d_ rj_ 

dz 8vr 



[-Bl + (x 2 -l)Bl\-p =0 



(33) 



and (32) by 



1 



[Bl-(x 2 -l)B 2 x ]+p = p eq . 



(34) 



As the equator is approached p steadily increases from the uniform value 
(implicit in the force- free assumption), which could in principle be zero, to 
the value p eq . Thus the model necessarily includes a thin domain with a 
non-force-free electromagnetic field. 

The limiting case, with the inequality sign in (30) replaced by equality, 
appears consistent with zero equatorial pressure p eq . However, condition 
(21) assumes that all non-electromagnetic forces, including inertial forces, 
are small compared with the dominant terms in the Lorentz force. I From 
(27), the rotation velocity is given by 



In a steady state, (11) yields the energy integral in a pressure-free system 
(MRWW, MS, Contopoulos 1995) 



showing that V would become infinite if = c/x, U p = c^J (x 2 — Tj/x; and 
from (35) this is equivalent to 



Thus the seemingly exceptional case, with zero equatorial pressure and so 
with the force-free equation holding all the way to the equator, in fact re- 
quires that the flow outside the equatorial zone be so highly relativistic that 
the neglected inertial terms are not small, so violating an essential condition 
for the force-free approximation to hold. We conclude that in all cases, the 
equilibrium conditions for this model will require at least a local breakdown 
in force-free conditions. In analogous non-relativistic problems, Lynden-Bell 
(1996) has pointed out that a thermal pressure is again required to balance 




(35) 



r[l — x(U^/c)] = constant, 



(36) 




(37) 
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magnetic pinching forces, exerted locally by an otherwise force-free field. 
However we again emphasize the important difference, that whereas in a 
non-relativistic problem, the electric stresses are normally smaller than the 
magnetic by the factor (U/c) 2 , in the present problem, beyond the 1-c the 
opposing electric stresses exceed the pinching poloidal field stresses by the 
factor x 2 , and equilibrium is possible only through the pressure exerted by 
the toroidal field. 

In the domain near the equator, with \B X /B Z \ 3> 1, the relativistic 
flow likewise has \U X /U Z \ S> 1, so that U 2 + C/| m c 2 , whence from (35), 
(remembering that B^ is negative), 

X-b^/(l + b 2 -X^) bx + y/(l+V-X*) 

U+/c= ^-g^ , U x /c= , (38) 

where b = \B ( j ) /B x \. (The algebraically allowed choice of the opposite signs 
before the two radicals would yield U^/c = 1 at x = 1, implying infinite T, 
and so is rejected). The outflow of angular momentum from the star in both 
hemispheres across a closed surface £ with local outward unit normal n is 
(e.g. Mestel 1999) 

r rPc 
- j (w^/47r)B p • ndS = -(c/a) 3 B? c J^ S(P)dP, (39) 

on use of (1), (5) and (22). 



4 The S(P) relation 

As noted in Section 2, if there do exist magnetospheric models that are 
everywhere non-dissipative, with the field force-free outside singular regions 
such as the equatorial sheet, then the relevant relation S(P) should emerge 
as part of the solution. An early model of the whole magnetosphere by 
Michel (1973a, 1991) has no dead zone, but a poloidal field that is radial all 
the way from the star to infinity. In our notation, the Michel field is 

P = Pc ( 1± (x 2 + z2) l/ 2 )' B * = ±P C(x 2 + z2) 3/2> B * = ±P c (x 2 + z2) 3/2> 

(40) 

with the negative sign applying to the northern hemisphere, and with the 
critical field line P = P c coinciding with the equator. From now on, just the 
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northern hemisphere is considered. Michel's S(P) relation is 



S 



p2 

-2P + — 



= ~P 



c 



(x 2 + z 2 )' 



X 



.2 



B<t> = ! 

x 



(41) 



so that 




Pc 



p 



) 



(x 2 + z 2 y/ 2 




{x 2 + z 2 y/ 2 ' 



(42) 



One can easily verify that (24) is satisfied. Note that the poloidal field (40) is 
radial and (away from the equator) independent of the spherical polar angle 
8, and so has no curl: equilibrium is maintained by balance between the 
electric force given by (16) and (17) and the force due to the toroidal field 
(41). However, with the sign change at the equator, there are again both 
toroidal and poloidal equatorial sheet currents that respectively maintain 
the jumps in B x and B^; also B 2 — {x 2 — l)B 2 = P 2 /x 4 when z = 0, so 
that by the above discussion, p eq > 0. Note also that dS/dP = on the 
critical line P c , but is negative for < P/P c < 1: in the Michel model, all 
the return current is in the equatorial sheet. 

In the generalized CKF picture, with the poloidal field forced to have a 
dipolar rather than a radial structure near the star, there is an associated 
dead zone within the 1-c, similar to that found in the MMWP, FM and 
MS fields, and analogous to that in the non-relativistic wind problem (e.g. 
Mestel and Spruit 1987). The dead zone terminates at the point (x c ,0); the 
value of x c (< 1) will be seen to be an extra parameter, fixing the global 
field structure. Within the dead zone the field lines close, crossing the 
equator normally so that the appropriate equatorial boundary condition is 
P z oc B x (x,0) = for x < x c . The dead zone is bounded by the separatrix 
field line P c . In the wind zone outside of P c , the wind flow is along the 
poloidal field, and so the equatorial boundary condition is P x oc B z = 0, 
but now extending through the 1-c from x c to oo, with again a charge- 
current sheet maintaining the discontinuous sign changes in E Z ,B X and B^. 
Thus the critical field line P(x, 0) = P c extending along the equator is the 
continuation of the separatrix between the wind and dead zones. 

Just outside the equatorial charge-current sheet and the separatrix, (24) 
holds all the way in from oo. At the intersection (1,0) with the 1-c and just 
outside the sheet, with P x oc B z = and no singularities in P xx and P zz , 
(24) requires that the constant value of SdS/dP along P c must be zero. But 
S(P C ) / 0, because as seen above, beyond the 1-c the pinching force exerted 
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by B,p is necessary for equilibrium; we therefore require, as in the Michel 
field, 

■dS" 



5pJ=° (43) 

on the critical field line P c . This condition is propagated inwards along P c on 
to its continuation, the separatrix P c between the wind and dead zones: the 
equatorial equilibrium conditions beyond the 1-c impose a constraint on the 
global S(P) relation. From (18), condition (43) requires that the poloidal 
volume current density j p falls to zero on P c , but the finite value for S(P C ) 
requires a poloidal current sheet at the equator beyond x c and along the 
separatrix within it. 



5 The domain within the light-cylinder 

We have emphasized that beyond the neutral point (x Cl 0), there is a finite 
thermal pressure on the equator, necessary for equilibrium. Likewise, it is 
not obvious that within the 1-c, the boundary condition on the field line 
separatrix between the wind zone (labelled 2) and the dead zone (labelled 
1) can be satisfied without a thermal pressure within the dead zone (see Fig. 

!)• 

It is in fact easy to generalize (21) and (24) to include a pressure gradi- 
ent: with gravity and inertia still negligible, 



47T I ZU 



.'aru\ 2 2B z /4tt\ 2 dS| „ 



(44) 

Thus p = p(P) - the constant pressure surfaces must coincide with the 
poloidal field lines. In normalized form, (44) becomes 

(1 X) dx* x dx +il X) dz*- V 2 X dp' (45) 

Prima facie, there is no obvious objection to the adoption of the simplest 
case, with p the same constant on all the field lines within the dead zone 
and zero in the wind zone, i.e. with a discontinuity in both p and S(P) on 
the separatrix. The equation for P within the wind zone then remains (24) 
(with signs reversed for convenience): 

(1 X W x dx +{1 X) dz*- V' (46) 
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Figure 1: Schematic magnetosphere 



and in the dead zone we have the same equation with S = 0. However, along 
the separatrix P c , extending inwards from the equatorial point (x c , 0), the 
equilibrium conditions require a discontinuity in B p as well as those in p and 
S{P). Writing 

m = -[(tx (B p /B p )]i (47) 
the unit normal to the separatrix, we require continuity of 

[- (8wp + E 2 + B 2 V + Bl)Sij + 2E i E j + 2B i B J -]n J - (48) 

with E given by (7). The components of (48) parallel to t and to the 
separatrix are automatically zero. The component normal to P c reduces to 
continuity of 

8Trp + B 2 [l-(aw/c) 2 ]+B 2 , (49) 

i.e. to 

8tt P + B 2 pl (l - x 2 ) = B 2 p2 (l - x 2 ) + B 2 2 . (50) 

It is convenient to normalize p in units of B 2 c /8ir. 

For x > x c , the wind zone 2 extends from the equator z = to z = oo; 
for x < x c , zone 2 extends from z c (x) defined by the separatrix 

P(x,z c ) = P c (51) 
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to z = oo. The separatrix function (51) is not known a priori but must 
emerge as part of the solution by iteration. We assume provisionally that 
at the point (x c ,0) (referred to as N), the poloidal field B p = both just 
outside and just inside the separatrix. From (50), the pressure of the toroidal 
field is balanced by the thermal pressure p, so that 

p = S 2 {P c )/x 2 c . (52) 

For x < x c , with use of (22) and (23), (50) then becomes 

(VPO 2 = (VP 2 ) 2 + S 2 (P c ) {1 ~f^ c) . (53) 

The discontinuity in |VP| grows from zero as x moves in from x c but will 
become a small fraction of |VP| for small x. 

When x c < 1, the function P has a simple analytical behaviour near 
the neutral point (x c , 0). The separatrix leaves (x c , 0) making an angle 
9 = 27r/3 with the outward-pointing equator. To leading order, in both the 
dead and wind zones, P xx + P zz = 0, which has the local solution 

P/P c = 1 + A 2 ,iR 3/2 sin(30/2), R 2 = z 2 + (x - x c ) 2 , (54) 

where the coefficients A21 apply respectively to the wind zone < 9 < 2tt/3 
and the dead zone 2-7r/3 < 9 < ir. It is seen that along 9 = 0, P = P c , and 
on 9 = ir, Pq cx P z = 0, as required. Across the separatrix z = V3(x c — x), 
P is continuous, while the jump condition (53) then yields 

A i = A * + 9P?x c (i-x2 c y (55) 

At this point, it is instructive to make a comparison with the analogous 
non-relativistic problem, in which the electric stresses are small by factors 
0(v/c) 2 and so are negligible, and also inertial forces are still neglected. 
Suppose that there is again a separatrix passing through the poloidal field 
neutral point at (x Cl 0). The balance equation across the separatrix is now 
continuity of 87rp + B 2 + B 2 ^ - the (1 — x 2 ) factors in (50) are replaced by 
unity. If again p is negligible in the wind zone, then the constant value of p 
along the separatrix within the dead zone is again fixed by the condition at 
the neutral point: p = S 2 (P c )/x 2 in normalized form. The balance condition 
is then S 2 (P c )(l - x 2 /x 2 ) + (VP 2 ) 2 = (VP!) 2 . The presence of the factor 
(1 — x 2 /x 2 ) in the S 2 term enables B p to be continuous (and zero) at the 
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neutral point (x c , 0), but with S 2 (P C ) still non-zero. (Clearly, in the non- 
relativistic problem, the numerical value of x c is of no significance.) 

By contrast, in the relativistic problem, near the 1-c the electric field 
strength approaches the poloidal magnetic field strength, so that the factor 
(1 — x 2 ) now appears multiplying both the (VP) 2 terms in (53). If now one 
were to take x c = 1, then the non-vanishing factor (1 — x 2 ) would cancel, 
and (53) would reduce to 

S 2 (P C ) + (VP 2 ) 2 = (VP 1 ) 2 . (56) 

along the separatrix. At the point (x c , 0), the simultaneous vanishing of VPi 
and VP2 would then require that S(P C ) = (and so also by (52) p = 0). 
But if beyond the 1-c the equatorial toroidal component Bs = —S(P c )/x 
were zero, then the balance condition could not be satisfied (cf. (32)). 

In fact, as discussed in the Appendix, the possible vanishing of S on 
the separatrix is appropriate to models with a radically different external 
equatorial boundary condition. In the present problem, the case x c = 1, 
which yields (56), may be incorporated provided that for this limiting case, 
we allow VPi / 0, i.e. the discontinuity in B p across the separatrix persists 
at the equator. If x c = 1 — e with e <C 1, then from (53), the equatorial B z 
can be zero on both sides of the separatrix, without requiring that S(P C ) = 0. 
However, at a neighbouring separatrix point x = x c — X with X < 1, the 
second term in (53) will have climbed from zero at N to S 2 {P C )[X/ (X + e)\ ~ 
S 2 {P C ) once X S> e. If solutions with S(P C ) / continue to exist as x c — > 1, 
then the separatrix balance condition will imply a steeper and steeper local 
gradient in B z , indeed tending to a discontinuity when x c = 1. 

For the limiting case x c = 1, an appropriate local model, valid near the 
critical point x = 1, z = 0, has |VP| 2 jumping by a constant value across 
the separatrix P = P j P c = 1 . New independent variables (r, t) and the 
dependent variable u are introduced, defined by 

x = l — rsint, z = rcost, r 2 = (x — l) 2 + z 2 ; P = 1 + u. (57) 

When r and u are small, (46) reduces to 

(1 - x 2 )(u xx + u zz ) - (1 + x 2 )u x /x = (58) 

in the dead zone. In the contiguous wind zone , the right-hand side is 
again SS', but as this will behave like u, the form (58) is appropriate on 
both sides of the separatrix. To leading order, this reduces to 

u„ + 2— + 1 \&aitut]t = 0. (59) 
r r z sin t 
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The appropriate boundary conditions are: 

(1) Ut = on t = 7r/2 - normal crossing of the equator; 

(2) u = on the separatrix, leaving the critical point at the as yet unknown 
angle t = to; 

(3) (u t /r) 2 jumps by a constant on t = tg. 

The solution of (59) satisfying these conditions is: 

Domain a < t < to : u = 0, 

Domain b to <t < ir/2 : u = ArQi[cost], (60) 

where Qi(X) is the Legendre function Qi = -\ + {X/2) log[(l+X)/(l-X)], 
and to, given by Qi(costo) = 0, is 33.5342°. The constant A then follows 
from u t /r = S(l). In Domain b, Q\[cost] = —1 — cost log(tant/2), and 
dQ\/d(cost) = cost/ sin 2 t — log(tant/2); whence at points near the 1-c (x 
close to unity), 

B x oc du/dz = — Alog(tant/2), B z oc — du/dx = ^4/sint, (61) 

yielding B x = 0, B z oc A = constant on the equator. 

An essential step in all the argument is the condition p = p(P), following 
from (44) . If instead the thermal pressure p were (illicitly) allowed to vary so 
as to balance the pressure B^ 2 = S 2 (P c )/x 2 exerted by the external toroidal 
field, then from (50) there would be no discontinuity in B p at the equator 
or indeed anywhere along the separatrix. 

It is appropriate also to re-emphasize just how crucial to the discussion 
is the outer domain equatorial boundary condition 

B z (x,0) =0 (62) 

for x > 1; for it is the consequent local equilibrium requirement S(P C ) > 
that implies a poloidal sheet current within the 1-c along the separatrix, 
leading to the conditions (50) and (53), with a non- vanishing B^- It is then 
clear that a non-zero p(P c ) is required for the balance condition to hold at 
the neutral point. 

6 Model construction 

The two basic parameters of a pulsar are clearly the angular velocity and 
the dipole moment of the neutron star. The spin-down rate of an active 
pulsar will depend on the strength of the circulating poloidal current. In 
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the Michel S(P) relation (41), the form with S ~ —2P valid near the poles 
derives from taking the current density flowing along the axis as equivalent 
to the G-J electron charge density moving with the speed c, a reasonable 
upper limit, and indeed one that must be closely approached if the electrons 
are to generate an electron-positron pair plasma near the star (cf. e.g. MS, 
Section 3). We consider only S(P) relations that behave like the Michel 
form for small P. 

The only other parameter introduced into the theory is the pressure in 
the dead zone, which we have taken as uniform. Prima facie, there should 
be a class of possible functions S(P), each fixing the global field function 
P and simultaneously the separatrix parameter P c and the limit x c of the 
dead zone, with the value of the required dead zone pressure p(P c ) given by 
(52). Equivalently, a choice of x c should have an associated class of functions 
S(P), each fixing the global P, the value of P c and of p(P c ), and the shape 
of the separatrix P(x,z) = P c between the star and x c . Intuitively, one 
expects an increase in p(P c ) to correspond to a decrease in x c . 

The functions P(x, z), S(P) are renormalized in terms of the value P c 
pertaining to the separatrix and its extension along the equator z = 0, 
extending from the point x c < 1, through the light-cylinder x = 1 to oo: 

P = P/P c , S = S/P c . (63) 

Although by convention the unnormalized P is negative, the new normal- 
ized form P is positive, and S is negative. (From now on the bar will be 
dropped, all quantities being assumed renormalized.) On the separatrix 
and its equatorial continuation, the normalized P is unity. The boundary 
condition at the origin then becomes P ~ (— 1/ P c )x 2 /(x 2 + z 2 ) 3 / 2 . 
The equation for P is 

(x 2 - 1)P XX + { -^±P x + (x 2 - 1)P„ = S%. (64) 
x dP 

Defining the spherical polar coordinates (r, 6) by z = rcosO, x = rsinO, 
it is to be expected that P ~ r n f{6) as r — ► oo, for some power n and 
function /. Such a homogeneous form is consistent with the equatorial 
condition P = 1 on 9 = ^tt, only if n = 0, giving a radial field. The 
appropriate boundary condition at infinity is thus dP/dr = 0. For large r, 
(64) takes the form x 2 V 2 P = SS'(P) which integrates once to require 

dP 

sine— = S(P) imposing S(0)=0. (65) 
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For the Michel field P = 1 — cos 6 and S = IP — P 2 , but a solution consistent 
with the dipole at the origin has a different function S(P). It follows that 
the radial field at infinity must differ from the Michel form. It is however 
found in the next section that the difference is not so great. 

6.1 Numerical Method 

The function P(x, z) satisfies an elliptic, quasilinear equation with a mixture 
of Neumann, Dirichlet and Robin boundary conditions. The PDE contains 
an unknown function S(P) which must be determined by a regularity con- 
dition across the singular line x = 1. There is a discontinuity in normal 
derivative across the curve P = 1 whose position is unknown a priori. The 
solution is driven by the dipole of strength 1/P C at the origin. A further 
numerical difficulty derives from the non-analytic behaviour of the solution 
near the point (x c ,0). 

The problem is determined by the two parameters P c and p, but for con- 
venience solutions are sought for fixed P c and x c . This reduces the amount 
by which the curve P = 1 requires adjustment during the solution. 

Equation (64) is discretised using a finite difference scheme on a rectan- 
gular grid (m5x, n5z) for 1 < m < M and < n < N, and P mn denotes the 
sought approximation to P at the grid point. The 1-c is at m = too, so that 
moSx = 1, and the computational domain is < x < x m , < z < z m where 
x m = Mbx and z m = N5z. 

The position of the dead-zone boundary P = 1 is marked by its intersec- 
tions with the grid lines. Away from this boundary and for to / too, second 
order centred difference approximations are used for P x , P xx and P zz , re- 
garding SS'(P) as a known source term. In fact, the dipole at the origin is 
subtracted out from P before differencing, and the dipole is differentiated 
analytically. 

Equations for P mo n and for P mn near the curve P = 1 are derived 
separately as discussed below. The resulting equations are written in a form 
appropriate for Jacobi iteration. Thus at each iteration new values of P mn 
are derived from the old values, for given S(P) and position of the curve 
P = 1. Once the new values of P mn are found, values for SS'(P) are found 
using 

SS (Pmon) — (-frno+ln Pmo — 1 n ) / Sx . 

These values are then interpolated using splines to give the new function 
SS'(P), which is integrated to give a new value of S(l). For values of P 
smaller than those reached along the finite length < z < z m the Michel 
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solution was used for S (P) . The intersection points of the curve P = 1 with 
the grid lines are then updated by requiring the jump in | VP| 2 to be correct. 
The internal and external values of |VP| are then used on the next iteration 
when approximating the derivatives at the grid points close to P = 1. The 
point (x ci 0) is kept fixed during the iteration, and a local balance there 
requires the system to converge to suitable values of p and S(l). 

Assuming a smooth crossing of the 1-c, P(x, z) is expanded near x = 1 
as a power series in (x — 1), which is substituted into (64). This gives the 
results 

SS'{P) = 2P X and AP XX + 2P ZZ = -{SS') X (66) 

on x = 1. This latter equation is discretised and used for m = m in place 
of (64). Finally, on the edge of the domain, the conditions P(0, z) = 0, 
P z (x, 0) = for x < x c while P(x, 0) = 1 for x > x c , and xP x + zP z = on 
x = x m and on z = z m are imposed. 

The above equations are solved iteratively. The numerical behaviour is 
improved using continuation techniques from a nearby solution, and by the 
use of under-relaxation. As well as the fundamental parameters P c and x c , 
the domain and z m and the grid size 5x and 5z affect the numerical 

solution. These latter values can be varied to ensure the computation is 
robust. 

No difficulties are encountered in the vicinity of the 1-c nor the dipole. 
However, when the 'front' P = 1 crosses one of the grid points in the course 
of the iteration some readjustment occurs and the convergence is slower. 
Similar effects occur in the numerical solution of Stefan problems with a 
liquid-solid interface. 

Solutions are found for some but not all parameter pairs (P c , x c ). For 
a fixed value of x c , the scheme converges only if P c lies in some interval 
Pi(x c ) < P c < P2(x c ). As x c — > 1, the values of P\ and P2 decrease, and 
become proportionally closer. If P c is too small the dead-zone becomes non- 
convex, bulging towards the 1-c away from the equator. It is believed that 
solutions cease to exist once P-lines attempt to cross the 1-c three times 
rather than once. If P c is too large, the curve P = 1 reaches the equator 
at a value of x < x c . Some dead-zone P-lines then cross the equator in 
the interval (x c , 1), and are not linked with the star and such solutions are 
deemed irrelevant. 

The left diagram in figure 2 portrays P(x, z) in the limiting case x c = 1 
with P c = 1.46, for which it is found S(l) — 0.919. An approximately linear 
variation with distance from (1, 0) within the dead zone is found, in keeping 
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Figure 2: Contours of P(x,z) Left: x c = 1, P c = 1.46; Middle: x c = 0.86, 
P c = 1.65; Right; x c = 0.45, P c = 2.9. 



with (60). With step-lengths 5x = 5z = 0.025, the value of P(l — 5x,0) 
predicts a value of A in (60) in agreement with the theoretical estimate to 
two significant figures. 

Solutions with x c < 1 are similar in appearance, although the local 
structure around the point (x Cl 0) differs as discussed above. It was not 
possible to distinguish the predicted separatrix angles 30° of (58.7) and 33° 
of (60). The centre of figure 2 shows the case x c = 0.86 P c = 1.65 for a larger 
domain, while the right depicts with a weaker dipole P c = 2.9 and 

smaller dead zone x c = 0.45, with equivalently an increase in the poloidal 
flux extending to infinity. 

The corresponding S(P) variations are plotted in figure 3, along with the 
quadratic Michel relation. For x c some distance from the 1-c S(P) is mono- 
tonic and the field on the 1-c diverges from the equator. As x c approaches 
unity, a region of negative S'(P) develops near P = 1, indicating that on 
the 1-c P x becomes negative so that the field is converging towards the equa- 
tor, although it becomes radial further out. The functional form of S(P) 
never seems to differ wildly from the Michel solution, but P(x, z) is affected 
substantially. This supports an observational point made at the outset of 
this paper, that magnetosphere currents may need very little alteration to 
achieve a different field configuration, in particular changes in the polar cap 
radius. Our solutions here correspond to polar caps which are respectively 
20% (x c = 1), 30% (x c = .86) and 70% (x c = .45) greater than the vacuum 
dipole case conventionally assumed in interpreting pulsar profile widths. 
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P 

Figure 3: S(P), for the cases in figure 2. Bottom to top: x c = 0.86, 1, 0.45, 
Michel. 

7 Discussion 

7.1 The equatorial boundary condition 

After considerable effort we have amended the pioneering study of CKF 
and constructed a class of models subject in the wind zone to the equatorial 
boundary condition B z = 0, i.e. with field lines in the wind zone never cross- 
ing the equator, as is indeed required by strict application of the 'perfect 
conductivity', dissipation- free condition. It has been emphasized that al- 
though the constructed magnetospheric field is force-free nearly everywhere, 
it is non-force-free in the the pinched equatorial zone where the field changes 
direction, and also along the separatrix between the wind and dead zones. 
Also, in general, as found earlier in CKF, the current returns to the star 
partly as a volume current, and partly as a sheet current along the equator 
and the separatrix. 

By definition, the global force-free field condition describes an electro- 
magnetically dominant system. In non-relativistic MHD, one would asso- 
ciate a poloidal field, drawn out so as to have a predominant x— component, 
rather with a kinetically-dominant system, with the Reynolds stresses at 
least locally stronger than the Maxwell stresses. Without such an out- 
ward pull, an equatorially pinching, quasi-radial field is likely to be unstable 
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against recurrent reconnection, spontaneously converting into a structure 
with field lines closing across the equator. 

Recall that in the analogous stellar problem (e.g. Mestel and Spruit 
1987, Mestel 1999, 2003), it is argued that the Alfvenic surface, on which 
the magnetic and kinetic energy surfaces are comparable, not only defines 
precisely the extent of 'effective corotation,' but is also a rough demarcation 
between the zones with respectively closed and open field lines. Suppose now 
that in studying the relativistic force-free equation, we consider models with 
a domain in which some field lines cross the 1-c and close beyond it. Gas in 
this domain cannot corotate with the star, but will flow into the equatorial 
zone. A hypothetical steady state will clearly require that the accumulating 
gas be able to flow outwards, crossing the field through a macro-resistivity. 

The most extreme example will have these field lines approaching the 
equator normally, so that the equatorial boundary condition B x = 0, holding 
within the 1-c, also replaces B z = 0, for some way beyond the 1-c. However, 
there is a further significant difference from the non-relativistic problem: 
a state with the equatorial boundary condition B x = 0, B z / appears 
to be electromagnetically unstable. If the field is perturbed so that there 
appear components ±B X on the ±z surfaces of the equatorial sheet, then 
the Maxwell stresses yield x-components (E X E Z + B x B z )/2tt; and provided 
coupling with the rotating magnetic star is maintained, this becomes (1 — 
x 2 )B x B z /2ir. Since B z < 0, when x > 1, the force acts to drive the frozen-in 
gas in the direction of B x , so increasing the imposed perturbation. The effect 
is again due to the electric part of the stress, which beyond the 1-c dominates 
over the poloidal magnetic field contribution. Thus there may exist steady 
states, alternative to those studied in this paper, in which the field at the 
equator is neither parallel nor perpendicular to the plane, but approaches 
it obliquely, and with the equatorial gas driven out by a combination of 
electrical and centrifugal force, and crossing the B z field component via a 
dynamically-driven macro-resistivity. 

Any such change in the equatorial boundary condition will react on the 
global solution of the force- free equation (cf. the Appendix). In MS (1994), 
two cases were studied, both subject to the special equatorial condition 
B x = 0: (a) the MMWP dead model, with S = 0, extrapolated beyond the 
1-c; and (b) the special live model, with S = —2P + 2P 2 /Pd, where Pj is the 
separatrix bordering the dead and wind zones. Note that in this tentative 
live model, S(Pd) = 0: all of the current returns as volume current, with no 
current sheet on the separatrix. In both these cases, it was found impossible 
to construct a solution that was smoothly continuous on the 1-c, and that 
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extended to infinity: the solutions became singular at x = 2 and x = 1.4 
respectively (cf Mestel 1999 and 2003, p. 547 and p. 612). 

For the second case, it was proposed that there exists near x = 1.4 
a thin volume domain (taken as cylindrical) in which E • B / 0. A pos- 
sible physical reason for the effective resistivity is the MHD instability of 
a local dominantly toroidal field (e.g. Begelman 1998). Prima facie, this 
allows construction of a solution that is well-behaved at both the 1-c and 
infinity. However, it is possible that the failure to find a global solution 
of the dissipation-free equation might turn out to be just a consequence of 
the assumed absence of a current sheet on the separatrix. This remains a 
problem for the future, along with a possible generalization to the cases with 
oblique approach of the field to the equator. A fully convincing treatment 
will require study of the gas and current flow in the equatorial zone. 

7.2 The inertial terms 

The solutions constructed according to the procedure of Section 6 may be 
regarded as appropriate generalizations of the Michel field, modified by im- 
position of the dipolar field on the star, with its associated dead zone. The 
parameters of the problem are assumed to yield the electromagnetic energy 
density dominating over the kinetic energy density, so that the deviation 
from force-free conditions is restricted to the thin equatorial domain. As 
the outflowing gas will be accelerated to high T-values, we need to check 
that the inertial terms remain small far beyond the 1-c. 

Refer back to the equations (10)-(14). We are interested now in the 
parameter domain in which the ideal MHD condition (13) remains an ade- 
quate approximation to (12), but the non-linear inertial term in (11) may 
no longer be negligible. We follow the standard treatment of a perfectly 
conducting, relativistic wind, flowing in the presence of the asymptotically 
radial Michel field and the associated toroidal field (Michel 1969; Goldreich 
and Julian 1970; Li and Melrose 1994; Mestel 1999, Section 7.9), using the 
notation of Melatos and Melrose. In addition to the kinematic relation (27), 
the continuity condition (10) has the steady state integral 

mnn = r/ (67) 

with m the electron/positron mass, and the equation of motion (11) yields 
the modified torque and energy integrals 

-wB^/Att + rjTmU^ = -0(P) /4vr (68) 
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and 

Tc 2 {l-awU (j> )=H{P). (69) 

(Recall that in this relativistic problem the electric force density p e E = 
— p e U x B/c is not negligible; however, it contributes to force balance across 
the field, but not to the integrals (68) and (69).) 

For simplicity, we consider just field lines near to the equator, with 

B m = <S>/zu 2 , (70) 

where $ is a measure of the poloidal flux crossing the 1-c. The principal 
results are the prediction of the asymptotic values 

rj ro ~c, CT"^ ~ 0, ~ -(anv/c)B w , (71) 

and 



a 2 <& 

47T7/C 3 



1/3 

= a 1 / 3 . (72) 



The force-free approximation will remain valid as long as the kinetic 
energy density is small compared with the electromagnetic, i.e. 

^ Tnmc 2 

ns wm <<A <73) 

Substitution from (67), (70), (71) and (72) yields 

n ~ a- 2/3 ~ r~ 2 . (74) 

The ±e-density n is written conveniently as the product of a pair-production 
multiplicity factor M and an estimate for the local Goldreich-Julian density: 
n = (M/e)(/9 e )Qj ~ M(aB p /2irec). Substitution of numbers yields 

r^[(10 7 /M)( J B sl2 ^ 6 /p2 )] i/3 (75) 

in terms of a stellar field B s = 10 12 G, stellar radius R s = 10 6 and pulsar 
rotation period Pi = lsec. Thus even with M = 10 3 , the predicted Tl < 1 
for the longest periods, so that the basic approximation leading to (21) 
and therefore the solutions of Section 6 - constructed subject to the crucial 
equatorial boundary condition B z = - appears to remain valid over the 
whole domain. 

The result is a self-consistency check for the present class of model, 
but also shows up a basic limitation to its applicability to a real pulsar; for 
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the condition (73) is - not surprisingly - also the condition that the energy 
flow at infinity be primarily by the Poynting flux rather than the kinetic 
flux. In fact, the wind flow results are highly model-dependent. Even with 
the same equatorial boundary condition, the predicted dominantly toroidal 
field at infinity may be dynamically unstable (e.g. Begelman 1998). Fur- 
ther, experience with other models with non-radial field lines beyond the 1-c 
(MRWW, FM, MS; Beskin et al 1993) suggests strongly that in an alterna- 
tive model with field lines approaching the equator obliquely (cf Sub-section 
7.1), particles streaming into the equator will have acquired much higher 
T- values. 

Finally, we emphasize that this paper has been concerned with the 
conditions for the construction of viable models. In practice, the particle 
pressure requirements in the equatorial and dead zones may very well put 
physical limitations on the model parameters - the extent x c of the dead 
zone and the associated separatrix P c . 
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Appendix: Fourier methods 
The present models 

Other studies of the pulsar magnetosphere have written the flux function 
P as a Fourier integral in z. In the illustrative MMWP model studied 
in MW and MP, in which the dead zone extends to the 1-c, P is written 
as a Fourier cosine integral, appropriate to a domain with the equatorial 
boundary condition B x oc dP/dz = 0. In the present paper, with the 
perfectly conducting wind domain extending from x c < 1 to oo and so with 
the equatorial boundary condition B z oc dP/dx = 0, the appropriate form 
for the normalized P-function in the domain x c < x < oo is the Fourier sine 
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integral 

2 r°° r°° 
P= — f(x, k) sin kzdk, f= Psinkzdz. (76) 
vr Jo jo 

This is a Fourier representation of the function P(x, z) in the north- 
ern hemisphere and —P(x,z) in the southern, with P(x,0±) = ±1. By 
Fourier's Theorem (e.g. Jeffreys and Jeffreys (1972)), at the discontinuity 
on the equator, the Fourier representation (76) should indeed have the value 
[(+1) + (— 1)]/2 = 0. However, this representation will exhibit the Gibbs 
phenomenon, and may make difficult the accurate construction of behaviour 
near the equator. 

It is convenient to introduce alsoF(x, k) defined by 

f=l+F(x,k). (77) 

Because of the boundary condition P = 1 on z = 0, the Fourier cosine 
transform of dP/dz is 

/ — coskzdz = -l+k / Psinkzdz = -l+kf(x,k) = kF(x,k). (78) 
Jo oz Jo 

Likewise, application of the Fourier sine transform to the force-free equation 
(64) yields 

(1 - x 2 )F" - (1 ^ 2) F / - (1 - x 2 )k 2 F = -g(x, k), (79) 
where as usual F' = F x , and 

roc 

g{x,k)= S—smkzdz. (80) 
Jo dP 

In general, the integral (80) can be converted into one over P: 

For a simple illustration, consider the Michel form (41) for S(P), for 
which (64) has the known solution (40) (again normalized with P c = 1), 
valid everywhere in the absence of a dead zone (i.e. with x c = 0). We quote 
the known Fourier transformation (AS, 9.6.25) 



JO 



COSkz T dz = ^ 2 ( —) V ^ {kx ) - . (82) 



{X 2 + Z 2)^+l/2 y 2 xj r^ + 1/2)' 
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Then substitution from (42) into (80) and from (40) into (76) and (77) yields 
after integration by parts 



F(x, k) = -xKi(b), g[x, k) = 2kx 2 K (kx), (83) 

whence (79) is seen to be satisfied from the known properties of Bessel 
functions. 

Now return to the problem with a dipolar field on the star and so with 
a finite dead zone, terminating at the point (x c ,0). It emerges from the 
computations reported in Section 6 that even with x c close to unity, the 
allowed S(P) relations do not differ much from the normalized Michel form 

S = -2P + P 2 . (84) 

The appearance of the neutral point (x c , 0) appears to be associated rather 
with the marked deviation from the radial structure for B p , found as the 1-c 
is approached from without. It is of interest to discuss this by application 
of the Fourier formalism. 

Having chosen a constructed model, and adopted the values found for 
P(l, z) and P x (l, z), one then moves inwards from x = 1, keeping P and P x 
continuous, but using the S(P) relation (84). One can write the solution of 
(79) as 

F(x,k) = b k F cf + Fpi. (85) 

where 6fc-F c / is a complementary function, satisfying (79) with zero r-h-s, and 
F p i is a particular integral, conveniently constructed to satisfy F p i(l, k) = 0. 
As in the MMWP problem, F c f is started off from x = 1 by the non-singular 
series 

F c/ = l + ^(l-x) 2 + ^(l-x) 3 + ...., (86) 

and is then continued inwards numerically. (Note that the function F c f for 
given k will remain unchanged during the whole operation.) 

With g(x, k) prescribed and F c f known, the particular integral is given 

by 



r-l 

Fpi = F c f I dx 



x U*[g(v)F cf {v)/v}dv) 



(87) 



(l + x)F* f (1-x) 
Near the 1-c, F p i is best computed from the series 
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found simply from (79) and its first derivative. 

With the particular choice for the P.I., F(l, k) = bhF c f(l, k)+F p i(l, k) = 
6^, so bk is fixed for given S(P) from the previously constructed solution 
for x > 1. For the present problem, from the given values of P on the 
1-c, continuity of F x and so of P x is ensured by the construction of F p i. 
Numerical continuation of F p i from (87) will depend on substitution into g 
of the values for P emerging from the inward integration. 

The model with x c = -86 was chosen for study. The inward integration 
is to be halted when a value x = x c is reached for which, using (78) together 
with (85), 



However, the numerical work yielded only partially satisfactory results. The 
preliminary integrations confirmed that the complementary function domi- 
nates - i.e., the solution is indeed determined primarily by the distribution 
of P on the 1-c - but the predicted value for .82. Attempts to 

bring x c closer to the value .86 found in the global integration of Section 6 
by including higher fc-values led to oscillating or unstable back transforms. 
We conclude that in a problem with the boundary condition P x (x, 0) = 0, 
requiring the Fourier sine integral (76), but with P(x, 0) = 1, the method is 
in general inappropriate for study of the field structure near the equator. 
The equatorial boundary condition B x = 

For the models of this paper, recall that in the inner domain, the dead 
zone terminates within the 1-c at the point N with coordinates (x c , 0) with 
x c < 1. Within the dead zone the field lines cross the equator normally, so 
that dP/dz = for x < x c , z = 0; whereas the condition that the wind zone 
be perfectly conducting enforces a fully open field crossing the 1-c, so that 
dP/dx = for x > x c , z = 0. This in turn yields a non-zero S both at the 
outer domain equator and on its continuation as the separatrix P c between 
the wind and dead zones, so that much of the current returns to the star as 
a sheet. Simultaneously, the conditions of equilibrium require a gas pressure 
both at the equator beyond x c and within the dead zone. 

In the somewhat more complicated model of MS there is again a force- 
free domain extending beyond the 1-c, but now the boundary condition 








(89) 



(90) 
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dP/dz oc B x (x,0) = is supposed to hold in the outer domain x > 1 
also. This is the appropriate approximation if e.g. a strong dynamically- 
driven macro-resistivity, allows the radial wind flow along the equator to 
cross the field lines. As in MW and MP, the solution for P can now be writ- 
ten as a Fourier cosine integral. The constraint (18) still holds away from the 
equator, but a local departure from torque-free conditions, associated with 
the flow of gas within the equatorial zone, allows S(P) to vary along the 
equator. Further, with B x (x,0) = 0, the arguments in Section 3 requiring 
a non-zero S at the equator and so along the separatrix between the wind 
and dead zones no longer hold. It is tentatively postulated that the current 
closure condition can now be satisfied by the simple choice S(P C ) = 0, al- 
lowing S(P) to go to zero continuously, without there needing to be a sheet 
current along P c . There is no pinched, high pressure equatorial sheet, nor 
is there a mandatory introduction of pressure into the dead zone. 

However, as noted in Section 7, these simplifications come at a price: 
for the case studied, with no current sheet on the separatrix, and with the 
equatorial boundary condition dP/dz oc B x = holding everywhere, a so- 
lution that is well-behaved and continuous at the 1-c blows up before it can 
reach infinity. In MP, it is shown that imposition of smooth continuity at 
the 1-c yields a singularity in the solution at less than two 1-c radii. Equiv- 
alently, one can apply standard asymptotic theory to the Fourier integral, 
demanding that the Fourier transform behave properly at infinity; there 
then appears an incompatibility of sign between solutions within and with- 
out the 1-c (Mestel 2001). Resolution of the dilemma may be as suggested 
in MP, with local breakdown in the simple plasma condition E • B = 0, not 
only on the equator, but also in a thin, dissipative volume domain, ideal- 
ized as a cylindrical shell symmetric about the rotation axis. Further work 
is needed to test whether instead, relaxation of the condition S{P C ) = 
will allow construction of a global solution of (64) subject to the condition 
B x (x,0)=0. 
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